function plotrocks(filenumber)
    N=readdata(0,1);
%     S=readdata(0,2);
%     T=readdata(0,3);
%     multigridlevels=readdata(0,4);
%     bound=readdata(0,5);
time=readdata(0,6);
%     dtnom=readdata(0,7);
dx=readdata(0,8);
dz=readdata(0,9);
mx=readdata(0,10);
mz=readdata(0,11);
%     mcs=readdata(0,12);
%     mct=readdata(0,13);
%     mtype=readdata(0,14);
%     mTemp=readdata(0,15);
%     mstressxx=readdata(0,16);
%     mstressxz=readdata(0,17);
%     mstrainxx=readdata(0,18);
%     mstrainxz=readdata(0,19);
%     mfinitestrain=readdata(0,20);
%     vx=readdata(0,21);
%     vz=readdata(0,22);
%     P=readdata(0,23);

mcs=zeros(1,N); mct=zeros(1,N);
xdivision=100; xplotcellsize=sum(dx)/xdivision;
x=0; for s=0:xdivision-1 ind=find(mx>=x & mx<x+xplotcellsize); mcs(ind)=s; x=x+xplotcellsize; end
zdivision=100; zplotcellsize=sum(dz)/zdivision;
z=0; for t=0:zdivision-1 ind=find(mz>=z & mz<z+zplotcellsize); mct(ind)=t; z=z+zplotcellsize; end
Ia=int32(find((-1).^mcs.*(-1).^mct>0));
Ib=int32(find((-1).^mcs.*(-1).^mct<=0));
clear('mcs','mct');

mtotalstrain=zeros(N,1); oldtime=time;

figure(2); clf;
for count=filenumber
    
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
    mtype=readdata(count,14);
    %mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
     %P=readdata(count,23);

    dt=time-oldtime; oldtime=time;
    x=zeros(1,S); z=zeros(1,T); 
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    figure(1);
    
    [grid,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,2,2);
    [Sgrid,Tgrid]=size(grid); surface=zeros(1,Sgrid); moho=zeros(1,Sgrid);
    surfaceind=zeros(1,Sgrid); mohoind=zeros(1,Sgrid); surface=zeros(1,Sgrid);
    for s=1:Sgrid
        depthslice=grid(s,:);
        ind=find(depthslice<1);
        [nothing,t]=max(zgrid(ind));
        surface(s)=zgrid(ind(t));
        surfaceind(s)=ind((t));
        ind=find(depthslice<2.99);
        [nothing,t]=max(zgrid(ind));
        moho(s)=zgrid(ind(t));
        mohoind(s)=ind(t);
        ind=find(depthslice<0.5);
        [nothing,indlow]=max(zgrid(ind)); indlow=ind(indlow);
        ind=find(depthslice>=0.5);
        [nothing,indhigh]=min(zgrid(ind)); indhigh=ind(indhigh);
        a=(grid(s,indhigh)-grid(s,indlow))/(zgrid(indhigh)-zgrid(indlow));
        b=grid(s,indlow)-a*zgrid(indlow);
        surface(s)=(0.5-b)/a;       
    end
    
    surface=surface-mean(surface(1:5));
    plot(xgrid,surface); set(gca,'ydir','rev'); axis([0 1000000 -3000 5000]);
    print(gcf, '-djpeg100',['..\outfiles\surface' num2str(count)]); clf;
    
    mtype(Ib)=mtype(Ib)+0.5;
    [grid,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,1,1);
    %[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; axis([0 1000000 0 400000]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); clear('Xgrid','Zgrid','grid','mtype');
    imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; axis([0 1000000 0 400000]); title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); clear('grid','mtype');
    print(gcf, '-djpeg100','-r600',['..\outfiles\rock' num2str(count)]); clf;
    
    mTemp=readdata(count,15);
    [grid,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
    %[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('Xgrid','Zgrid','mTemp');
    imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('mTemp');
    print(gcf, '-djpeg100','-r600',['..\outfiles\temp' num2str(count)]); clf;
    [Sgrid,Tgrid]=size(grid); q=zeros(1,Sgrid); mohotemp=zeros(1,Sgrid);
    for s=1:Sgrid %this loops assumes that the grid is the same used when finding the moho and surface
        t=surfaceind(s);
        q(s)=3*(grid(s,t+5)-grid(s,t+1))/(zgrid(5)-zgrid(1)); %assuming k=3
        t=mohoind(s);
        mohotemp(s)=grid(s,t);
    end
    clear('grid');
    [ax,nothing,nothing]=plotyy(xgrid,q,xgrid,mohotemp); set(ax(1),'xlim',[0 1500000]); set(ax(2),'xlim',[0 1500000]); set(ax(1),'ylim',[3e-2 25e-2],'ytick',3e-2:1e-2:25e-2); set(ax(2),'ylim',[0 800],'ytick',0:50:800);
    print(gcf, '-djpeg100','-r600',['..\outfiles\surfaceq-mohotemp' num2str(count)]); clf;

    mstrainxx=readdata(count,18); mstrainxz=readdata(count,19);
    mtotalstrain=mtotalstrain+dt*(mstrainxx.^2+mstrainxz.^2).^0.5;
    [grid,xgrid,zgrid]=bilingrid(mtotalstrain,mx,mz,dx,dz,1,1);
    %[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('Xgrid','Zgrid','grid');
    imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('grid');
    print(gcf, '-djpeg100','-r600',['..\outfiles\totalstrain' num2str(count)]); clf;
    
    [grid,xgrid,zgrid]=bilingrid((mstrainxx.^2+mstrainxz.^2).^0.5,mx,mz,dx,dz,1,1); clear('mstrainxx','mstrainxz');
    %[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('Xgrid','Zgrid','grid');
    imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('grid');
    print(gcf, '-djpeg100','-r600',['..\outfiles\strain' num2str(count)]); clf;

    mfinitestrain=readdata(count,20);
    [grid,xgrid,zgrid]=bilingrid(mfinitestrain,mx,mz,dx,dz,1,1); clear('mfinitestrain');
    %[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); pcolor(Xgrid,Zgrid,grid); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('Xgrid','Zgrid','grid');
    imagesc(xgrid,zgrid,grid'); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]); clear('grid');
    print(gcf, '-djpeg100','-r600',['..\outfiles\plasticstrain' num2str(count)]); clf;

% mtype=readdata(count,14); P=readdata(count,23); mstressxx=readdata(count,16); mstressxz=readdata(count,17); mTemp=readdata(count,15);
% zP=zeros(1,T)+0.5*dz(1); xP=zeros(1,S)+0.5*dx(1);
% for t=1:T-2 zP(t+1)=zP(t)+0.5*(dz(t)+dz(t+1)); end
% zP(T)=zP(T-1)+dz(T-1);
% for s=1:S-2 xP(s+1)=xP(s)+0.5*(dx(s)+dx(s+1)); end
% xP(S)=xP(S-1)+dx(S-1);
% [ZP,XP] = meshgrid(zP,xP); P=P(1:end-1,1:end-1); XP=XP(1:end-1,1:end-1); ZP=ZP(1:end-1,1:end-1);
% [gridP,xgrid,zgrid]=bilingrid(P(:),XP(:),ZP(:),dx,dz,1,1); clear ('XP','ZP','P');
% [grid,xgrid,zgrid]=bilingrid((mstressxx.^2+mstressxz.^2).^0.5,mx,mz,dx,dz,1,1); clear('mstressxx','mstressxz');
% [gridtype,xgrid,zgrid]=bilingrid(mtype,mx,mz,dx,dz,1,1); clear('mtype'); [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
% gridtype=round(gridtype); dens=zeros(size(gridtype));
% ind=find(gridtype<1); dens(ind)=1; ind=find(gridtype>=1); dens(ind)=2800; ind=find(gridtype>=2); dens(ind)=3000; ind=find(gridtype>=3); dens(ind)=3300;
% [gridTemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1); clear('mTemp');
% dens=dens-3.25e-5.*gridTemp; clear('gridTemp');
% hydrostat=zeros(size(gridtype));
% for s=1:S
%     for t=2:T
%         hydrostat(s,t)=hydrostat(s,t-1)+9.82*0.5*(dens(s,t)+dens(s,t-1))*dz(t-1);
%     end
% end
% grid=grid-gridP+0*hydrostat; %sigmaxx=sigmaxxdev-(P-hydrastat)
% intgstress=zeros(1,S);
% for s=1:S
%     intgstress(s)=trapz(zgrid,grid(s,:));
% %     plot(zgrid,grid(s,:));
% %     plot(zgrid,gridP(s,:),'r');
% %     plot(zgrid,hydrostat(s,:),'g');
% %     plot(zgrid,gridP(s,:)-hydrostat(s,:),'c');
% %     plot(zgrid,grid(s,:)-gridP(s,:)+hydrostat(s,:),'k');
% end
% clear('grid','gridP','hydrostat','dens','gridtype');
% %clf; plot(intgstress);
% %print(gcf, '-djpeg100',['..\outfiles\stressdepth' num2str(count)]);
% avintgstress=trapz(xgrid,intgstress)/xgrid(end);
% figure(2); hold on; plot(time/(60*60*24*365*1e6),avintgstress,'r.'); figure(1);
% %pcolor(Xgrid,Zgrid,gridP-hydrostat); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 400000]);

    
%     figure(1); plot(mx(I0(I0a)),mz(I0(I0a)),'y.','markersize',1); hold on; plot(mx(I0(I0b)),mz(I0(I0b)),'g.','markersize',1);
%     plot(mx(I1(I1a)),mz(I1(I1a)),'b.','markersize',1);
%     plot(mx(I1(I1b)),mz(I1(I1b)),'c.','markersize',1);
%     set(gca,'ydir','rev'); hold off; axis equal; title([num2str(time/(60*60*24*365*1e6)) ' Ma ']); axis([0 1000000 0 250000]);
%     print(gcf, '-djpeg100',['..\outfiles\a' num2str(count)]);
% 
%     [gridstrain,xgrid,zgrid]=gridmarkers(mfinitestrain,mx,mz,mcs,mct,x,z,S*2,T*2); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(2); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     %[gridstrain,xgrid,zgrid]=gridmarkers1(mfinitestrain,mx,mz,mcs,mct,x,z); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(2); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     print(gcf, '-djpeg100',['..\outfiles\plasticstrain' num2str(count)]); 
%     clear ('gridstrain','xgrid','zgrid'); clf;
%     I=find(mtype==0); mstrainxx(I)=0; mstrainxz(I)=0; %resetting air strain
%     [gridstrain,xgrid,zgrid]=gridmarkers((mstrainxx.^2+mstrainxz.^2).^0.5,mx,mz,mcs,mct,x,z,S*2,T*2); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(5); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     %[gridstrain,xgrid,zgrid]=gridmarkers1((mstrainxx.^2+mstrainxz.^2).^0.5,mx,mz,mcs,mct,x,z); [Xgrid,Zgrid] = meshgrid(xgrid,zgrid); figure(5); pcolor(Xgrid,Zgrid,gridstrain); shading flat; set(gca,'ydir','rev'); axis equal; colorbar; axis([0 1000000 0 250000]);
%     print(gcf, '-djpeg100',['..\outfiles\' num2str(count)]);
%     clear ('gridstrain','xgrid','zgrid');
 
end
end